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Abstract 

We simulate a finite system of N confined electrons with inclusion of the 
Darwin magnetic interaction in two- and three-dimensions. The lowest energy 
states are located using the steepest descent quenching adapted for velocity 
dependent potentials. Below a critical density the ground state is a static 
Wigner lattice. For supercritical density the ground state has a non-zero 
kinetic energy. The critical density decreases with N for exponential confine- 
ment but not for harmonic confinement. The lowest energy state also depends 
on the confinement and dimension: an antiferromagnetic cluster forms for 

harmonic confinement in two dimensions. 
PACS numbers: 05.70.Fh, 41.20.-q, 64.70.-p, 36.40.-c 
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I. INTRODUCTION 



The use of fast modern computers has made it increasingly easy to investigate many- 
dimensional finite systems and to study dynamical quantities like time to equipartition, sym- 
metry break in finite systems and quantities of thermodynamic interest The numerical 
studies have also revealed unexpected features of finite systems like cluster formation |^-|^ 
and lack of equipartition In this paper we study numerically the energetic quenching 
of a classical electron gas with inclusion of the magnetic interaction The main motiva- 
tion to study this system is the physics it describes and unfortunately the dynamics of this 
finite system is not easy to study because the equations of motion are algebraic-differential 
in character [|l^]. Because of this we explored numerically only the quenching motion, which 
involves integration of simple ordinary differential equations. Like the Coulomb interaction, 
the magnetic interaction is a long-range interaction. The long range nature of the Coulomb 
interaction has been widely explored in molecular dynamics and there is a large literature 
of numerical simulations and various special techniques were invented to deal with the long 



range nature of the interaction |TT|. On the other hand, the magnetic interaction has been 
much less studied. 

The magnetic interaction appears naturally in classical Maxwell electrodynamics: The 
lowest-order retardation and magnetic effects (or order (f/c)^) can be described in terms of 
electron variables only as a velocity-dependent interaction. This approximation was orig- 
inally proposed by Darwin [|l^] to obtain a lagrangian which bears his name. The Darwin 



lagrangian is much-used in atomic physics |]T3|,|14| where it is known as Darwin-Breit inter- 
action in its quantised form. The Darwin lagrangian includes the lowest order correction 
to the electric field of a moving charge and the lowest order magnetic field (the Biot-Savart 
term). Apart from its traditional domain of atomic physics the Darwin lagrangian has been 
used to model slightly relativistic plasmas |[T5Hl^ and even models of superconductivity 
and stellar magnetic fields 0. These ideas extend the range of Darwin lagrangian from 
conventional relativistic correction: the Darwin corrections are also important in the low- 



2 



energy (nonrelativistic) regime. One reason for that is to be found in the long-rangedness 
of the Darwin force and the fact that it is badly screened, unhke the Debye screening of 



the Coulomb force [|1^]. These two factors compensate for the weakness of the Darwin in- 
teraction at low velocities and make it a player in certain circumstances The Darwin 
lagrangian has also been used as an unfolding of the scale invariant degenerate Coulom- 
bian interaction to estimate long-time-scale dynamical effects in atomic physics [T^. Last, 
the Darwin lagrangian represents the first correction to the Coulomb interaction lagrangian 
in a series expansion of the Tetrode-Fokker lagrangian of the Wheeler- Feynman action at 
a distance theory. This relativistic lagrangian theory has been shown to have interesting 



magnetic consequences pO| that deserve further study. 

Experimentally, aggregates of electrons with the same sign of charge can be confined 
for long times using suitably chosen static electric and magnetic fields and form so-called 
nonneutral plasmas (For a recent review see l21|]). Differently from neutral plasmas (i.e. 
composed of electrons with both sign of charge in approximately equal numbers), the non- 
neutral plasmas can attain thermal equilibrium and cooled to low enough temperatures to 
form liquid and crystal-like states. These plasmas under various confinement geometries 
have been extensively studied, particularly in the last decade. Generally the interelectron 
interaction is taken to be Coulombic or Yukawa. 

The study of the possible low-energy states of a system may shed some light on the 
possibility of a phase transition. We chose to study the magnetic long-range effects in a 
confined electron gas and in the neighborhood of its lowest energy state. We confine a 
system of charged electrons in two and three dimensions by use of a background field, taken 
either to be of a harmonic or exponential form. The electrons interact by Coulomb plus 
Darwin forces. For the purely Coulombian repulsion, it is known that the lowest energy 
state of this system is a Wigner lattice (triangular-like in two dimensions), which we also 
find with our quenching techniques p2| , p3| . We show that the long range effects depend on 
N and on a single parameter f3, which is usually of the order of 10^^ — 10^'^ for attainable 
physical densities. We find that low-A^ systems need an artificially larger value of f3 for the 



long-range effects to be important. The assumption is that the critical value goes down with 
A^, and because of the practical impossibility of simulating systems with millions of electrons, 
we investigate finite systems for an artificially higher value of (3 and extrapolate the scaling 
properties to the large-A^ case. This assumption holds for the exponential confinement but 
not for the harmonic confinement, where the critical (3 does not go down with A^. 

This paper is divided as follows. The next section introduces the model and defines 
the quenching procedure. The numerical results for the 2-D systems are described in the 
section 3; two subsections correspond to types of confining potential. The section 4 deals 
with the 3-D systems, again divided into subsections. The papers ends with a discussion in 
section 5. 



II. DARWIN LAGRANGIAN AND THE NATURAL QUENCHING 

We consider N electrons in two or three dimensions interacting via the Coulomb re- 
pulsion plus the velocity-dependent Darwin magnetic interaction |]9|,p!2[| and confined by a 



one-electron potential Vc{t) of the positive background. The lagrangian for this system can 
be written as 



TV TV 1 TV 



where and Vi are the position and the velocity of the ith electron and r'jj = fi — fj, 
Tij = \ fij\, iij = Tij/rij is the unit vector pointing from the ith to the jth electron, e is the 
electronic charge, and c is the velocity of light. The first term is the kinetic energy of 
the system and the second term is the Coulomb energy. The next term is the Darwin Vd 
which can be expressed in terms of the vector potential A as 

Vn = -^Tv,.A, (2) 
zc • 



with 



A,(r,t;)^eE^^|^^^- (3) 



As the lagrangian is time-independent, there is an associated energy constant which evaluates 
to [0] 

E = Y: + Y^va^{r, ^) + E ^ + E yc{rd. (4) 

i i<j *J i 

Notice that this energy is not of the minimal coupling type 

E=J:l{p^~Ar+v, (5) 

which is only the case when the magnetic field is external to the system, and consequently the 
vector potential is velocity independent. In the present situation, because of the internal 
fields, the state of minimal energy is not always the zero velocity case anymore, as the 
conditions for the Bohr- Van Laufen theorem are not satisfied P, p!B| , P^ (This theorem states 
that a velocity independent vector potential does not affect the partition function if the 
energy is of the minimal coupling form (^)). 

The form of the equations can be simplified by using scaled units: a length scale, given 
by the average interelectron separation R, scales positions as x Rx (in this units the gas 
has a density equal to one). Time is scaled as dt — > uodr where Uq = e'^/mR^. In these 
units the energy scales as — > m{uJoRyE with 

^ = E^^t + E;^+/5'E^^^^%^^^^^^ + E^c(rO. (6) 

i i<j U i<j *J i 

The parameter jS"^ in the above equation is defined as 

- 5, (7) 

where = e'^/mc^ is the classical electronic radius and the interelectron distance R in 
2-D is given by i? = l/\/n and R = n~^^^ in 3-D. For some real physical situations: The 
conduction band in metals forms a 3-D degenerate plasma with typical densities of ~ 10^^ 
cm~^, which gives for the value of ~ 10^^. The highest density physical plasma is 
found in the interior of white dwarf stars, corresponding to an electron density of n ~ 10^^ 



cm ^ which gives /5 ~ 10 ^ ||2^ 



Once the systems studied here are rotationally invariant, Noether's theorem determines 



a constant of motion p9[ for them 



i=l i<j '^'^ij 

with Cjj = fij/rij, as before. For 2-D this constant is a vector perpendicular to the plane 
and for 3-D the 0(3) symmetry determines a constant vector by the above formula. This 
constant can be interpreted simply as the sum of the mechanical angular momenta plus the 
field angular momentum. 

To look for the minimum energy states of the velocity dependent A^-body system, we 
adapt a numerical procedure analogous to the steepest descent quenching using what we 
name the natural quenching vector field. We check that for potential systems this procedure 
produces the static crystalline arrangement of electrons known as the Wigner lattice p2| , p3| . 
We now define the natural quenching vector field, which is constructed from the differential 
of the expression for the energy constant given by equation (P). We start from a random 
initial condition and integrate it as a function of a "quenching parameter" by the following 
gradient equations: 

dfi dE dvi dE , . 

ds dfi ' ds dvi 

It is easy to see that along this gradient motion the energy always decreases, as the parameter 
derivative of the energy evaluates to 

^ = _ VI— 12- VI— 12 dni 

ds ^'afj ^^dvi^ ■ ^ ^ 

Along with the numerical quenching from random initial conditions it is necessary to 
use the relativistic form of the kinetic energy. Otherwise, we observe that some electrons 
acquire an enormous kinetic energy during quenching, creating an enormous nonphysical 
internal field that still decreases the total energy. Of course, for such large velocities the 
Darwin approximation breaks down and the whole lagrangian describes nonphysical effects. 



as discussed in reference [^. In all our numerical experiments we check that the electron 



energies were never relativistic in the final quenched state, which guarantees that the Darwin 
approximation is vahd. Last, to gain some understanding of how the above quenching 
procedure can find states with nonzero velocity, let us examine equation (|^) for the velocities, 
which read as 

dvi^ _ _-j _ ^2 ^ t/j + eij{vj.eij) ^^^^ 
ds 2rij 

Notice that on the right hand side we have a linear function of the velocities, defining a 
linear matrix M(rj,/?^). For (3'^ = 0, this matrix is minus the identity and the velocities are 
all quenched down to zero. Above a critical value of jS"^, this matrix can have negative and 
zero eigenvalues, and it is not possible to quench the velocities to zero anymore, which is 
the cause of the nonzero velocity states we find. The critical point (3"^ can also be located 
by an alternative analytical method: Consider equation (|ll|) for the velocity-quenching. For 
= 0, the eigenvalues of M are all degenerate and equal to one. Taking the electron 
coordinates to be those of the static Wigner lattice (which can be obtained for (3"^ = 0), 
one can diagonalize M numerically and find all its eigenvalues. The critical is that for 
which the minimum eigenvalue of M crosses zero i.e. the minimum eigenvalue of M is just 
negative. It can be seen then that in this case the quenching will decrease the energy while 
increasing the velocities along the negative eigenvector directions. We have diagonalized M 
in the neighborhood of the Wigner lattice and it is satisfactory that the critical calculated 
by the matrix method agrees with the values obtained by quenching. 

III. NUMERICAL RESULTS FOR CIRCULAR DISK GEOMETRY 

A. Harmonic confinement 

We consider first a system of electrons in 2-D, confined by the field of a uniformly 
charged circular disk of positive charges, of radius Rd scaled units, and with the electronic 
density of one electron per squared scaled unit (A^ = vri?^). For this system the potential 
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of the uniformly charged disk of positive background can be calculated analytically p8[ and 
for r < Rd it is approximated by 

3/2 

Mr) = -2i{nN) + (12) 

We have explicitly included the negative constant to properly account for the electrostatic 
interaction with the positive background. One still needs to add the self-energy of the 
positive background ^/ttN^^'^/S to expression @ in order to get the total electrostatic energy. 
We seek to determine the lowest-energy states of this system by employing the natural 
quenching technique described above, and we integrate equations @ numerically with an 
6/7 Runge-Kutta embedded integrator pair. By quenching from different initial conditions 
we can hope to obtain insight into the character of the ground state. The natural quenching 
is performed for the disk system for various values of the parameter The electrons are 
started from a triangular lattice, distorted slightly in a random manner, with velocities 
uniformly (and randomly) distributed upto a certain maximum value. A square lattice type 
initial configuration is also used and the same final result is found. The system is quenched 
until a steady state appears to have been reached. To check if we actually attain a global 
minimum state and not merely a local minimum, we slightly heat the obtained configuration 
and quench it again. By these means we are confident that our ground states are at least 
quahtatively correct. 

These simulations are done for 225 to 1600 electrons in the disk. In all cases it is observed 
that below a certain value of the parameter jS"^ the ground state is the static Wigner lattice, 
independent of the But above the critical a new type of ground state is obtained. 
This state has nonzero kinetic energy with a striking nonuniform distribution of velocities. 
The electrons with large velocities are confined to an antiferromagnetic cluster in the center 
of the disk. The configuration in the position space remains visibly triangular-like. The 
electrons in the central cluster have velocities aligned in a manner to minimize the Darwin 
energy (Fig.|T]). This parallel and antiparallel orientation of the velocities succeeds in lowering 
the energy of the nonstatic configuration below the static Wigner lattice. 
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The critical parameter value (31 decreases with the increasing A^, the number of electrons 
in the disk, but it is rather a weak dependence. For the 225 electrons disk, (31 ~ 0.72 
while it is approximately 0.71 for the 1600 electrons disk. The relative size of the cluster 
slowly decreases with increasing A^: we quantify it as following. A electron i belongs to the 
cluster if (3'^vf > 0.01. The quantity is generally about 0.06 for the fastest electron. 
At (3"^ = 0.75, this criterion yields the fraction of the cluster electrons to be 0.31 for the 225 
electron system and monotonously decreasing to 0.26 for A^ = 900 and to 0.23 for A^ = 1600. 
However it appears that the cluster remains equally hot, independent of A^, i.e., the average 
kinetic energy per cluster electron does not depend on A^ at constant (3^. Futher increasing 
the value of beyond Pi causes rapid increase both in the size and the temperature of 
the cluster. The ground-state energy continues to decrease as is increased beyond Pi 



The quenching runs do not always yield the ground state. Often, in particular for larger 
sized disks, an imperfectly aligned higher energy state is obtained. The local order is of 
the same type as the true ground state but on a larger scale two or more regions of the 
local order have an mismatch, analogous to grain-boundaries in a polycrystalline material 



The effect of choice of the confining potential may be studied by considering a different 
potential. To this end, we replace the harmonic potential by an exponential potential. 



Here, as before Rd = ^{N/tt) and r^y = 0.5 in scaled units. The quenchings are performed 
in a manner identical to that described above, starting from a randomly distorted triangular 
lattice. The critical (3l is obtained, separating the static and non-static ground-states. The 
values of Pi are rather lower and they decrease appreciably with increasing A^. We find 




(Fig. §). 



B. Exponential confinement 



Vc = Voexp((r - Rd)/rw)- 
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131 ^ 0.52 for = 225 to about 0.42 for = 1600. These (and the intermediate = 484 
and = 900) values may be fitted to a power-law (31 ~ A^" with a ~ —0.11. 

The character of the ground-states is also affected. The edge (i.e. the surface) of the 
static lattice is no longer triangular-like but is composed of two ring-like layers. Above 
the velocity distribution is highly inhomogeneous: the kinetic energy is concentrated in the 
two edge layers (Fig. |). Increasing beyond (3l does not result in more electrons acquiring 
kinetic energy but merely increases kinetic energy of the edge electrons. 

IV. NUMERICAL RESULTS FOR SPHERICAL GEOMETRY 

A. Harmonic confinement 

In this section we consider A^ electrons confined by the field of a homogeneous positively 
charged sphere of radius Rd- The potential (for r < R^) can be calculated exactly, 

„9 27r r, 
Vc = -2^Rl + —r\ 

which follows immediately from Gauss law of electrostatics. To obtain the total electrostatic 
energy we must again add the self-energy contribution of the background ^NRj to the 
expression (P). As before, the electron density is taken to be one electron per scaled unit. 
Quenchings are performed for A^=216-1000 electron systems. Simulations are started from a 
randomly distorted cubic lattice, while the velocities are initialized in the manner previously 
described for the disk geometry. As in the disk geometry, a critical separates static 
ground states from the nonstatic ground states. The critical is slightly smaller in three 
dimensions — it varies from ^ 0.67 for A^ = 216 to about 0.66 for A^ = 1000. 

But the character of the ground states is very different from that obtained for the disk 
geometry. The electrons are arranged in a multiple ring-like structure around the center of 
the ball. These rings possess sharp boundaries and their number grows with A^. For example, 
the 216 electrons system has 4 such rings (including the cluster of central electrons) while 7 
rings are visible for A^ = 1000 system. 
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This ring-like structure persists beyond j]^. Though, the velocity distribution is not 
homogenous, a distinct cluster of hot electrons is not formed. The velocities project radially 
outwards (Fig. ^) . The electrons in the small central cluster have smaller than average kinetic 
energy. The kinetic energy is fairly shared between the other rings but the distribution 
between electrons in a particular ring is non-uniform. 

The different ordering in two and three dimensions has striking effects: in 3-D the Dar- 
win interactions between neighboring electrons are highly repulsive and the lowering of the 
energy is provided by the Darwin interactions of distant electrons. Whereas, in the 2-D disk 
geometry, the Darwin interactions between nearest neighbor electrons lower the Darwin 
energy. 

B. Exponential confinement 

In this subsection we describe the results of the simulation of 3-D electron gas confined 
by an exponential potential 

Vc = Voe^^{{r-Rd)lrw). 

Here Rd is defined by the relation = ^-R^, and r^/ = 1.5. The static and nonstatic ground 
states are again obtained, below and above respectively. Though, in contrast with the 
harmonic confinement, the (3^ decreases with increasing with the power-law (3^ ~ A^~*^'^^. 
The values of f3^ range from 0.74 for = 216 to 0.50 for N = 1000. 

The lattice has a ring-like structure again but with some differences. First, the number 
of rings is smaller: = 1000 system has only 4 rings compared with 7 with harmonic 
potential. Second, the central cluster is absent. The distribution of kinetic energy is also 
different from the harmonic case. The kinetic energy per electron increases as one goes 
outwards and for A^ = 1000 the inner two rings have almost no kinetic energy. 
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V. DISCUSSIONS AND CONCLUSIONS 



It would be natural to integrate the supercritical symmetry-breaking states of low-energy 
and nonzero angular momentum to see the time dependence of the angular momentum [Q, 
but this is not an easy numerical job. The existence of zero and negative eigenvalues of 
M in (|TT]) signal the onset of complex dynamical behavior for this system: The lagrangian 
equations of motion that follow from (|T]) are 

Q. _)_ ^2 ^ + eijjo-j.eij) _ dVc _ y-v Cij 
+/^^ ^[(^i ■ + (l^il^ - ■ Cij - 2vi ■ Vj)ei,j\. (13) 



Notice that this is an algebraic-differential equation of order zero |T0|, and the linear matrix 



sitting on the left side is the same exact matrix M that appeared in (|ri|). The numerical 
procedure to integrate this equation is delicate: If the matrix has maximal rank, which is 



the case for low values of /5, the integrator RADAU 30] can be used. If the matrix has 



only one zero eigenvalue, then the integrator DASSL |[TD| can be used, as long as the matrix 



does not lose rank along the trajectory, which is the case for rare initial conditions only. 
The general case where the matrix's rank is lesser than 2N — 1, or even worse if it loses 
rank along the trajectory, then one is faced with a rich system which could generate a very 
complex dynamics. As a matter of fact, in 1976, one of the earliest studies of this many-body 
system declared it intractable [|15| and a coarse grained field approximation was developed 



to study it, which later became a plasma simulation technique [|T7|. We could integrate 
initial conditions with a very low (3, and there we found that the angular momentum is an 
approximate constant, as it should be for the /5 = case, according to equation (||). For 
intermediate values of (3 we were able to perform the numerical integration of the dynamics 
using DASSL, and we find that the total angular momentum does not change sign for very 
long time scales. In the supercritical situation, where it would be of interest to study the 
dynamics, we find that the time steps of RADAU quickly go to zero due to the criticality 
of the matrix. An integration method has still to be developed to simulate this interesting 
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dynamics. 

In three dimensions, because of Gauss law, charge neutrahty causes the total electrostatic 
energy to be an extensive quantity. As a matter of fact, our data show that E/N tends 
to a constant value in 3-D. On the other hand our data for 2-D indicates that the total 
electrostatic energy is non-extensive. 

At constant the Darwin lowering of energy is smaller for larger N systems in 2-D. 
This can also be inferred from the observation made in section 3A that the fraction of 
hot electrons goes down with increasing A^. Thus our results indicate that as limA^ — >■ 
oo at constant electron density the Darwin lowering vanishes. This conclusion holds for 
harmonically confined systems for which (31 is almost independent of A^ and hence different 
sized systems can be compared at constant 

The magnetic lowering does not decrease with large A^ in 3-D. This makes the question 
of a proper thermodynamic limit more subtle. In a certain sense, the Darwin interaction 
merely renormalizes the electronic charge. But this would interfere with the cancellation 
of the background charge and hence jeopardize the thermodynamic limit |^. We feel that 
more numerical and analytical work is needed to resolve the question of the thermodynamic 
limit for the Darwin lagrangian. 

Acknowledgements: We acknowledge discussions with A. Castelo, F. C. Alcaraz and J. 
P. Rino. 
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FIGURES 

FIG. 1. The ground-state of harmonically confined disk with N = 900 at = 0.75. The 
electrons are located at the tails of the arrows; the lenght of the arrows is proportional to the 
magnitude of the electron velocities obtained after quenching. The direction of an arrow gives the 
angle of the corresponding velocity vector. Scaled coordinates. 

FIG. 2. (a) Ground-state energy of harmonically confined disk with N = 225 vs. /3^. The 
critical is slightly less than 0.72. (b) Ground-state energy of harmonically confined electron gas 
in 3-D with N = 216 vs. The critical /^^ is slightly less than 0.67. 

FIG. 3. A higher-energy minimum of harmonically confined disk with N = 900 at = 0.75. 
The arrows are made in the manner described in the previous figure. 

FIG. 4. The ground-state of exponentially confined electron gas in 2-D at = 0.45. The 
arrows are drawn as in Fig. 1. 

FIG. 5. The ground-state for harmonically confined system in three dimensions with N = 216 
at (3"^ = 0.74. The velocities project radially outwards. Arrows are made in the similar manner to 
the two-dimensional case. 
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